#!/usr/bin/env python

#  Original Algorithm/Concept:
#   Francisco Gambon
#     Servicio Gallego de Salud, Unidad de Inmunologia
#     Hospital de Meixoeiro, 36210, Vigo, Spain
#
#  Update Implementation:
#   David Olivieri
#     E.S. Exeñeria Informatica
#     Universidad de Vigo, 32004, Ourense, Spain.



import sys
import getopt, sys
#import numpy as np
from operator import itemgetter, attrgetter
import time
import re


from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
from Bio import Motif
from Bio.SeqRecord import SeqRecord
from Bio import SeqFeature
from Bio.Graphics import GenomeDiagram


class VregionsExtract:
    def __init__(self, inFile=None, outVfile=None,
                 outVRfile=None, geneFormat=None):

        print "VregionsExtract: inside Vregion constructor"
        text="inside constructor"
        strand=1
        self.inFile=inFile
        self.outVfile=outVfile
        self.outVRfile=outVRfile
        self.geneFormat=geneFormat
        self.add_MotifSeqs()
        self.ofile = open(self.outVfile,"w")


    def __call__(self):
        #print "inside call"
        self.add_MotifSeqs()
        self.start()


    def runVR(self):
        self.add_MotifSeqs()
        self.start()

    def add_MotifSeqs(self):
        text="inside add_Motif"
        self.m = Motif.Motif(alphabet=IUPAC.unambiguous_dna)
        self.m.add_instance(Seq("YYC",self.m.alphabet))
        self.m.add_instance(Seq("YFC",self.m.alphabet))
        self.m.add_instance(Seq("YLC",self.m.alphabet))
        self.m.add_instance(Seq("YIC",self.m.alphabet))
        self.m.add_instance(Seq("YHC",self.m.alphabet))
        self.m.add_instance(Seq("TFC",self.m.alphabet))

    def get_VregionFB2(self, Lexon):

        def run_fastExtrV(Record,iNum_Sequence,strand):

            if strand == 1:
                Sequence=Record.seq
                StringR = "  "
            else:
                Sequence=Record.seq.reverse_complement()
                StringR = "R "

            Label=Record.description
            pa_ant = -200
            ID = 1

            termino = Sequence.find('CAC',Lexon)
            while termino != -1:
                inicio = termino-Lexon
                vx = Sequence[inicio:termino]
                qcola=[]
                for t in re.finditer("AG", str(vx)[0:40]):
                    qcola.append(t.start())
                qcnt=0
                for q in qcola:
                    qcnt+=1
                    s = q + 4
                    vs = vx[s:-1]

                    if 'TGG' in vs[99:135] :
                        p = vs.translate(to_stop=True)
                        pa = termino-3*len(p)-4
                        if len(p) > 85:
                            z = p[33:45]
                            z2 = p[18:30]
                            z3 = p[-12:-1]
                            if 'C' in z2 and 'W' in z:
                                for pos,seq in self.m.search_instances(z3):
                                    pos>1
                                    if pa > pa_ant + 200:
                                        pa_ant = pa
                                        pa=((1-strand)/2)*len(Record)+strand*pa     
                                        recordV=SeqRecord(p, id = str(strand) + "|" + str(pa) + "|"
                                                          + str(ID), description = "V" + StringR + "region :: " + Label)

                                        self.ofile.write(recordV.format("fasta"))
                                        ID += 1
                                        break  


                termino = Sequence.find("CAC", termino+2)
                if termino > len(Record): break



        t0=time.time()
        print 'forward'
        irec = 0
        for record in SeqIO.parse(self.inFile, "fasta"):
            irec += 1
            run_fastExtrV(record,irec,1)

        print 'reverse(POSICIONES CORREGIDAS)'
        irec = 0
        for record in SeqIO.parse(self.inFile, "fasta"):
            irec += 1
            run_fastExtrV(record,irec,-1)


        dt=time.time()-t0
        print "FINAL TIME=",dt
        self.ofile.close()  

## ---------------MAIN ----------------------------------
if __name__ == '__main__':

    ## Allow for command line args and fixed.
    try:
        #V =  VregionsExtract(sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4])
        outVRfile='outVRregions'
        geneFormat='fasta'
        V =  VregionsExtract(sys.argv[1], sys.argv[2], outVRfile, geneFormat)
    except IndexError:
        inFile='gorila.txt'
        outVfile='outVregions'
        outVRfile='outVRregions'
        geneFormat='fasta'
        V = VregionsExtract(inFile, outVfile, outVRfile, geneFormat)

    #V.get_VregionFB()

    Lexon=320
    V.get_VregionFB2(Lexon)






